week 2: linear model and causal inference

geocentric models

Workspace setup:

library(here)
library(tidyverse)
library(cowplot)
library(brms)
library(tidybayes)
library(patchwork)

Model “recipes”

  1. Recognize a set of variables to work with. (Data and parameters.)
  2. Define each variable either in terms of the other variables OR in terms of a probability distribution.
  3. The combination of variables and their probability distributions defines a joint generative model that can be used to simulate hypothetical observations and analyze real ones.

Here’s an example:

\[\begin{align*} y_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \beta x_i \\ \beta &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \\ x_i &\sim \text{Normal}(0,1) \\ \end{align*}\]

Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

  • \(W\) is the observed count of water.
  • \(N\) is the total number of tosses.
  • \(p\) is the proportion of water on the globe.

The whole model can be read as:

The count \(W\) is distributed binomially with sample size \(N\) and probability \(p\). The prior for \(p\) is assumed to be uniform between 0 and 1.

Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

Estimating the posterior using brms

Last week, we used grid approximation to estimate the posterior distribution. In the video you watched for today, McElreath moves on to something called QUADRATIC APPROXIMATION. It’s good to understand what that’s doing, but you and I are moving right along to MCMC. We won’t go into the details of what’s happening for a few weeks, but let’s start with the code to estimate the model.

Let’s say we tossed the globe 9 times and observed 6 waters:

m1 <-
  brm(data = list(w = 6),
      family = binomial(link = "identity"),
      w | trials(9) ~ 0 + Intercept,
      prior(beta(1, 1), class = b, lb = 0, ub = 1),
      iter = 5000, warmup = 1000, seed = 3, chains=1,
      file = here("files/models/m21.1"))
1
Data can be a data frame or a list. Just make sure the variable names match what’s in your formula.
2
How you assume your outcome variable is distributed.
3
The formula for your outcome.
4
Priors for every parameter in your model. Here, we only have one parameter, so we only need 1 prior. This happens to be a flat prior between 0 and 1.
5
Some choices about how we want our model to run. We’ll go more into this later.
6
These models can take a long time to run. You can automatically save the output to a file; when you do this, the next time you run this code, it won’t actually estimate the model, but will instead pull the output from your stated file. Be WARNED: if you change the data or the model code, it will NOT restimate your model until you delete the file.

Sampling from the posterior

Grid approximation gave us the calculated probability of each possible value of our parameter, \(p\). But our method of conducting bayes will no longer give us such a neat solution. Here’s how you get the posterior distribution for \(p\):

samples_from_post = as_draws_df(m1)
samples_from_post
# A draws_df: 4000 iterations, 1 chains, and 3 variables
   b_Intercept lprior lp__
1         0.66      0 -2.8
2         0.66      0 -2.8
3         0.56      0 -2.9
4         0.36      0 -4.5
5         0.68      0 -2.8
6         0.62      0 -2.8
7         0.44      0 -3.7
8         0.67      0 -2.8
9         0.83      0 -4.0
10        0.75      0 -3.1
# ... with 3990 more draws
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
samples_from_post %>%  
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water")

Posterior predictive distribution

ppd = posterior_predict(m1)
ppd
        [,1]
   [1,]    6
   [2,]    9
   [3,]    7
   [4,]    3
   [5,]    6
   [6,]    6
   [7,]    5
   [8,]    6
   [9,]    7
  [10,]    5
  [11,]    9
  [12,]    6
  [13,]    9
  [14,]    7
  [15,]    6
  [16,]    7
  [17,]    4
  [18,]    5
  [19,]    5
  [20,]    9
  [21,]    0
  [22,]    7
  [23,]    3
  [24,]    8
  [25,]    7
  [26,]    6
  [27,]    6
  [28,]    6
  [29,]    8
  [30,]    7
  [31,]    8
  [32,]    8
  [33,]    5
  [34,]    3
  [35,]    7
  [36,]    3
  [37,]    5
  [38,]    7
  [39,]    7
  [40,]    9
  [41,]    7
  [42,]    5
  [43,]    6
  [44,]    3
  [45,]    5
  [46,]    4
  [47,]    4
  [48,]    6
  [49,]    5
  [50,]    6
  [51,]    7
  [52,]    9
  [53,]    7
  [54,]    4
  [55,]    7
  [56,]    8
  [57,]    5
  [58,]    3
  [59,]    5
  [60,]    3
  [61,]    3
  [62,]    4
  [63,]    3
  [64,]    4
  [65,]    4
  [66,]    4
  [67,]    4
  [68,]    6
  [69,]    5
  [70,]    5
  [71,]    7
  [72,]    5
  [73,]    7
  [74,]    6
  [75,]    7
  [76,]    6
  [77,]    6
  [78,]    7
  [79,]    7
  [80,]    4
  [81,]    5
  [82,]    9
  [83,]    7
  [84,]    2
  [85,]    5
  [86,]    5
  [87,]    4
  [88,]    6
  [89,]    6
  [90,]    6
  [91,]    6
  [92,]    5
  [93,]    3
  [94,]    4
  [95,]    4
  [96,]    7
  [97,]    7
  [98,]    5
  [99,]    8
 [100,]    8
 [101,]    5
 [102,]    5
 [103,]    5
 [104,]    6
 [105,]    8
 [106,]    6
 [107,]    3
 [108,]    7
 [109,]    6
 [110,]    4
 [111,]    6
 [112,]    3
 [113,]    7
 [114,]    5
 [115,]    5
 [116,]    7
 [117,]    3
 [118,]    6
 [119,]    6
 [120,]    7
 [121,]    6
 [122,]    7
 [123,]    6
 [124,]    6
 [125,]    6
 [126,]    7
 [127,]    8
 [128,]    2
 [129,]    4
 [130,]    3
 [131,]    4
 [132,]    1
 [133,]    3
 [134,]    3
 [135,]    7
 [136,]    7
 [137,]    2
 [138,]    6
 [139,]    5
 [140,]    4
 [141,]    7
 [142,]    7
 [143,]    8
 [144,]    5
 [145,]    5
 [146,]    9
 [147,]    7
 [148,]    7
 [149,]    6
 [150,]    5
 [151,]    5
 [152,]    7
 [153,]    4
 [154,]    6
 [155,]    7
 [156,]    5
 [157,]    7
 [158,]    6
 [159,]    8
 [160,]    7
 [161,]    8
 [162,]    7
 [163,]    6
 [164,]    5
 [165,]    5
 [166,]    5
 [167,]    4
 [168,]    3
 [169,]    5
 [170,]    7
 [171,]    7
 [172,]    7
 [173,]    8
 [174,]    8
 [175,]    8
 [176,]    6
 [177,]    4
 [178,]    6
 [179,]    7
 [180,]    6
 [181,]    4
 [182,]    4
 [183,]    5
 [184,]    7
 [185,]    5
 [186,]    6
 [187,]    7
 [188,]    9
 [189,]    7
 [190,]    4
 [191,]    7
 [192,]    5
 [193,]    3
 [194,]    2
 [195,]    5
 [196,]    2
 [197,]    5
 [198,]    5
 [199,]    6
 [200,]    7
 [201,]    4
 [202,]    5
 [203,]    3
 [204,]    7
 [205,]    9
 [206,]    8
 [207,]    8
 [208,]    8
 [209,]    8
 [210,]    8
 [211,]    9
 [212,]    7
 [213,]    8
 [214,]    2
 [215,]    5
 [216,]    7
 [217,]    6
 [218,]    7
 [219,]    4
 [220,]    6
 [221,]    5
 [222,]    7
 [223,]    3
 [224,]    9
 [225,]    5
 [226,]    7
 [227,]    6
 [228,]    6
 [229,]    6
 [230,]    6
 [231,]    3
 [232,]    4
 [233,]    5
 [234,]    4
 [235,]    6
 [236,]    4
 [237,]    7
 [238,]    4
 [239,]    4
 [240,]    5
 [241,]    7
 [242,]    6
 [243,]    5
 [244,]    4
 [245,]    6
 [246,]    7
 [247,]    9
 [248,]    5
 [249,]    8
 [250,]    7
 [251,]    7
 [252,]    6
 [253,]    7
 [254,]    8
 [255,]    8
 [256,]    6
 [257,]    7
 [258,]    7
 [259,]    5
 [260,]    7
 [261,]    4
 [262,]    3
 [263,]    7
 [264,]    6
 [265,]    6
 [266,]    7
 [267,]    4
 [268,]    4
 [269,]    5
 [270,]    7
 [271,]    3
 [272,]    3
 [273,]    3
 [274,]    6
 [275,]    4
 [276,]    4
 [277,]    6
 [278,]    9
 [279,]    7
 [280,]    7
 [281,]    6
 [282,]    5
 [283,]    9
 [284,]    3
 [285,]    6
 [286,]    7
 [287,]    6
 [288,]    5
 [289,]    6
 [290,]    7
 [291,]    6
 [292,]    5
 [293,]    3
 [294,]    6
 [295,]    8
 [296,]    7
 [297,]    5
 [298,]    6
 [299,]    7
 [300,]    7
 [301,]    7
 [302,]    5
 [303,]    2
 [304,]    5
 [305,]    0
 [306,]    3
 [307,]    3
 [308,]    3
 [309,]    5
 [310,]    3
 [311,]    3
 [312,]    6
 [313,]    7
 [314,]    4
 [315,]    5
 [316,]    6
 [317,]    6
 [318,]    3
 [319,]    1
 [320,]    7
 [321,]    8
 [322,]    6
 [323,]    8
 [324,]    7
 [325,]    8
 [326,]    8
 [327,]    7
 [328,]    6
 [329,]    4
 [330,]    6
 [331,]    6
 [332,]    7
 [333,]    8
 [334,]    3
 [335,]    7
 [336,]    3
 [337,]    5
 [338,]    8
 [339,]    5
 [340,]    6
 [341,]    4
 [342,]    7
 [343,]    4
 [344,]    5
 [345,]    8
 [346,]    6
 [347,]    7
 [348,]    2
 [349,]    5
 [350,]    8
 [351,]    7
 [352,]    5
 [353,]    6
 [354,]    8
 [355,]    5
 [356,]    8
 [357,]    4
 [358,]    8
 [359,]    6
 [360,]    6
 [361,]    6
 [362,]    5
 [363,]    7
 [364,]    7
 [365,]    6
 [366,]    7
 [367,]    6
 [368,]    6
 [369,]    8
 [370,]    7
 [371,]    6
 [372,]    9
 [373,]    7
 [374,]    4
 [375,]    8
 [376,]    8
 [377,]    6
 [378,]    4
 [379,]    4
 [380,]    4
 [381,]    6
 [382,]    4
 [383,]    4
 [384,]    5
 [385,]    4
 [386,]    7
 [387,]    6
 [388,]    5
 [389,]    7
 [390,]    4
 [391,]    6
 [392,]    8
 [393,]    6
 [394,]    8
 [395,]    6
 [396,]    2
 [397,]    8
 [398,]    5
 [399,]    5
 [400,]    8
 [401,]    3
 [402,]    8
 [403,]    7
 [404,]    5
 [405,]    9
 [406,]    4
 [407,]    3
 [408,]    6
 [409,]    7
 [410,]    5
 [411,]    6
 [412,]    8
 [413,]    6
 [414,]    4
 [415,]    5
 [416,]    6
 [417,]    5
 [418,]    4
 [419,]    2
 [420,]    6
 [421,]    8
 [422,]    6
 [423,]    4
 [424,]    6
 [425,]    5
 [426,]    4
 [427,]    3
 [428,]    6
 [429,]    4
 [430,]    5
 [431,]    5
 [432,]    6
 [433,]    8
 [434,]    6
 [435,]    6
 [436,]    4
 [437,]    5
 [438,]    4
 [439,]    3
 [440,]    5
 [441,]    5
 [442,]    5
 [443,]    5
 [444,]    5
 [445,]    5
 [446,]    8
 [447,]    5
 [448,]    5
 [449,]    3
 [450,]    4
 [451,]    1
 [452,]    3
 [453,]    2
 [454,]    4
 [455,]    5
 [456,]    8
 [457,]    5
 [458,]    8
 [459,]    5
 [460,]    7
 [461,]    7
 [462,]    6
 [463,]    9
 [464,]    4
 [465,]    8
 [466,]    6
 [467,]    3
 [468,]    3
 [469,]    5
 [470,]    7
 [471,]    8
 [472,]    7
 [473,]    6
 [474,]    7
 [475,]    7
 [476,]    4
 [477,]    4
 [478,]    9
 [479,]    6
 [480,]    7
 [481,]    8
 [482,]    5
 [483,]    7
 [484,]    8
 [485,]    9
 [486,]    6
 [487,]    7
 [488,]    5
 [489,]    2
 [490,]    6
 [491,]    6
 [492,]    4
 [493,]    6
 [494,]    5
 [495,]    7
 [496,]    5
 [497,]    8
 [498,]    7
 [499,]    6
 [500,]    6
 [501,]    4
 [502,]    8
 [503,]    5
 [504,]    3
 [505,]    4
 [506,]    7
 [507,]    7
 [508,]    3
 [509,]    5
 [510,]    4
 [511,]    2
 [512,]    9
 [513,]    7
 [514,]    7
 [515,]    5
 [516,]    8
 [517,]    3
 [518,]    8
 [519,]    6
 [520,]    9
 [521,]    3
 [522,]    6
 [523,]    5
 [524,]    7
 [525,]    6
 [526,]    7
 [527,]    9
 [528,]    6
 [529,]    4
 [530,]    2
 [531,]    3
 [532,]    4
 [533,]    8
 [534,]    8
 [535,]    9
 [536,]    4
 [537,]    6
 [538,]    8
 [539,]    8
 [540,]    2
 [541,]    7
 [542,]    7
 [543,]    7
 [544,]    5
 [545,]    7
 [546,]    6
 [547,]    7
 [548,]    5
 [549,]    4
 [550,]    5
 [551,]    5
 [552,]    6
 [553,]    7
 [554,]    6
 [555,]    8
 [556,]    6
 [557,]    8
 [558,]    9
 [559,]    8
 [560,]    8
 [561,]    4
 [562,]    4
 [563,]    5
 [564,]    3
 [565,]    4
 [566,]    8
 [567,]    8
 [568,]    8
 [569,]    4
 [570,]    5
 [571,]    7
 [572,]    7
 [573,]    6
 [574,]    5
 [575,]    5
 [576,]    5
 [577,]    3
 [578,]    3
 [579,]    4
 [580,]    8
 [581,]    6
 [582,]    9
 [583,]    6
 [584,]    6
 [585,]    5
 [586,]    7
 [587,]    4
 [588,]    5
 [589,]    5
 [590,]    5
 [591,]    5
 [592,]    5
 [593,]    7
 [594,]    6
 [595,]    8
 [596,]    5
 [597,]    4
 [598,]    5
 [599,]    6
 [600,]    2
 [601,]    6
 [602,]    6
 [603,]    6
 [604,]    7
 [605,]    5
 [606,]    5
 [607,]    6
 [608,]    8
 [609,]    3
 [610,]    4
 [611,]    4
 [612,]    7
 [613,]    5
 [614,]    2
 [615,]    6
 [616,]    5
 [617,]    7
 [618,]    8
 [619,]    8
 [620,]    7
 [621,]    7
 [622,]    5
 [623,]    6
 [624,]    5
 [625,]    3
 [626,]    4
 [627,]    1
 [628,]    4
 [629,]    4
 [630,]    8
 [631,]    8
 [632,]    4
 [633,]    7
 [634,]    5
 [635,]    6
 [636,]    8
 [637,]    7
 [638,]    4
 [639,]    7
 [640,]    7
 [641,]    2
 [642,]    5
 [643,]    6
 [644,]    8
 [645,]    6
 [646,]    7
 [647,]    8
 [648,]    6
 [649,]    6
 [650,]    8
 [651,]    6
 [652,]    9
 [653,]    8
 [654,]    5
 [655,]    2
 [656,]    6
 [657,]    4
 [658,]    1
 [659,]    7
 [660,]    7
 [661,]    5
 [662,]    5
 [663,]    4
 [664,]    8
 [665,]    6
 [666,]    7
 [667,]    7
 [668,]    4
 [669,]    6
 [670,]    8
 [671,]    7
 [672,]    7
 [673,]    6
 [674,]    7
 [675,]    6
 [676,]    8
 [677,]    6
 [678,]    5
 [679,]    4
 [680,]    3
 [681,]    7
 [682,]    8
 [683,]    4
 [684,]    4
 [685,]    3
 [686,]    5
 [687,]    5
 [688,]    8
 [689,]    7
 [690,]    5
 [691,]    4
 [692,]    6
 [693,]    4
 [694,]    2
 [695,]    6
 [696,]    5
 [697,]    2
 [698,]    6
 [699,]    7
 [700,]    6
 [701,]    8
 [702,]    3
 [703,]    8
 [704,]    8
 [705,]    0
 [706,]    3
 [707,]    8
 [708,]    7
 [709,]    6
 [710,]    3
 [711,]    5
 [712,]    1
 [713,]    7
 [714,]    4
 [715,]    8
 [716,]    7
 [717,]    8
 [718,]    6
 [719,]    5
 [720,]    8
 [721,]    6
 [722,]    5
 [723,]    7
 [724,]    9
 [725,]    6
 [726,]    9
 [727,]    8
 [728,]    7
 [729,]    8
 [730,]    2
 [731,]    8
 [732,]    7
 [733,]    4
 [734,]    8
 [735,]    4
 [736,]    4
 [737,]    7
 [738,]    5
 [739,]    7
 [740,]    7
 [741,]    6
 [742,]    4
 [743,]    5
 [744,]    4
 [745,]    6
 [746,]    7
 [747,]    6
 [748,]    6
 [749,]    7
 [750,]    6
 [751,]    6
 [752,]    3
 [753,]    7
 [754,]    6
 [755,]    7
 [756,]    7
 [757,]    8
 [758,]    7
 [759,]    9
 [760,]    9
 [761,]    7
 [762,]    7
 [763,]    6
 [764,]    2
 [765,]    6
 [766,]    6
 [767,]    6
 [768,]    8
 [769,]    3
 [770,]    7
 [771,]    6
 [772,]    7
 [773,]    7
 [774,]    5
 [775,]    5
 [776,]    8
 [777,]    4
 [778,]    3
 [779,]    3
 [780,]    5
 [781,]    5
 [782,]    6
 [783,]    6
 [784,]    7
 [785,]    4
 [786,]    8
 [787,]    3
 [788,]    7
 [789,]    6
 [790,]    5
 [791,]    6
 [792,]    8
 [793,]    7
 [794,]    5
 [795,]    7
 [796,]    6
 [797,]    9
 [798,]    8
 [799,]    7
 [800,]    7
 [801,]    7
 [802,]    8
 [803,]    6
 [804,]    5
 [805,]    8
 [806,]    9
 [807,]    8
 [808,]    8
 [809,]    3
 [810,]    9
 [811,]    7
 [812,]    7
 [813,]    7
 [814,]    4
 [815,]    5
 [816,]    6
 [817,]    6
 [818,]    7
 [819,]    7
 [820,]    8
 [821,]    3
 [822,]    7
 [823,]    7
 [824,]    6
 [825,]    7
 [826,]    3
 [827,]    6
 [828,]    6
 [829,]    4
 [830,]    3
 [831,]    4
 [832,]    4
 [833,]    2
 [834,]    5
 [835,]    8
 [836,]    4
 [837,]    3
 [838,]    4
 [839,]    5
 [840,]    2
 [841,]    2
 [842,]    8
 [843,]    7
 [844,]    4
 [845,]    6
 [846,]    6
 [847,]    8
 [848,]    4
 [849,]    7
 [850,]    6
 [851,]    6
 [852,]    5
 [853,]    9
 [854,]    6
 [855,]    1
 [856,]    4
 [857,]    7
 [858,]    7
 [859,]    7
 [860,]    8
 [861,]    6
 [862,]    7
 [863,]    7
 [864,]    6
 [865,]    6
 [866,]    8
 [867,]    4
 [868,]    8
 [869,]    6
 [870,]    4
 [871,]    6
 [872,]    3
 [873,]    8
 [874,]    8
 [875,]    6
 [876,]    5
 [877,]    7
 [878,]    5
 [879,]    4
 [880,]    6
 [881,]    6
 [882,]    4
 [883,]    4
 [884,]    5
 [885,]    6
 [886,]    6
 [887,]    4
 [888,]    7
 [889,]    5
 [890,]    6
 [891,]    6
 [892,]    8
 [893,]    6
 [894,]    6
 [895,]    1
 [896,]    1
 [897,]    2
 [898,]    6
 [899,]    8
 [900,]    5
 [901,]    4
 [902,]    5
 [903,]    1
 [904,]    2
 [905,]    5
 [906,]    4
 [907,]    3
 [908,]    3
 [909,]    5
 [910,]    5
 [911,]    2
 [912,]    6
 [913,]    7
 [914,]    7
 [915,]    6
 [916,]    6
 [917,]    4
 [918,]    4
 [919,]    8
 [920,]    7
 [921,]    5
 [922,]    2
 [923,]    4
 [924,]    4
 [925,]    7
 [926,]    5
 [927,]    6
 [928,]    7
 [929,]    7
 [930,]    5
 [931,]    4
 [932,]    6
 [933,]    5
 [934,]    6
 [935,]    5
 [936,]    7
 [937,]    4
 [938,]    4
 [939,]    4
 [940,]    6
 [941,]    4
 [942,]    4
 [943,]    1
 [944,]    1
 [945,]    2
 [946,]    6
 [947,]    7
 [948,]    6
 [949,]    7
 [950,]    6
 [951,]    6
 [952,]    8
 [953,]    3
 [954,]    5
 [955,]    3
 [956,]    2
 [957,]    4
 [958,]    6
 [959,]    9
 [960,]    6
 [961,]    4
 [962,]    7
 [963,]    3
 [964,]    9
 [965,]    5
 [966,]    3
 [967,]    5
 [968,]    2
 [969,]    8
 [970,]    9
 [971,]    5
 [972,]    5
 [973,]    4
 [974,]    3
 [975,]    7
 [976,]    7
 [977,]    9
 [978,]    8
 [979,]    8
 [980,]    7
 [981,]    7
 [982,]    3
 [983,]    5
 [984,]    5
 [985,]    8
 [986,]    6
 [987,]    3
 [988,]    4
 [989,]    5
 [990,]    2
 [991,]    6
 [992,]    7
 [993,]    9
 [994,]    6
 [995,]    7
 [996,]    8
 [997,]    2
 [998,]    3
 [999,]    2
[1000,]    8
[1001,]    5
[1002,]    5
[1003,]    5
[1004,]    4
[1005,]    4
[1006,]    6
[1007,]    4
[1008,]    4
[1009,]    9
[1010,]    8
[1011,]    2
[1012,]    4
[1013,]    8
[1014,]    2
[1015,]    7
[1016,]    8
[1017,]    7
[1018,]    8
[1019,]    4
[1020,]    6
[1021,]    8
[1022,]    7
[1023,]    9
[1024,]    8
[1025,]    7
[1026,]    4
[1027,]    6
[1028,]    6
[1029,]    6
[1030,]    4
[1031,]    8
[1032,]    6
[1033,]    8
[1034,]    3
[1035,]    5
[1036,]    8
[1037,]    7
[1038,]    8
[1039,]    7
[1040,]    8
[1041,]    9
[1042,]    5
[1043,]    8
[1044,]    5
[1045,]    5
[1046,]    7
[1047,]    7
[1048,]    6
[1049,]    6
[1050,]    7
[1051,]    2
[1052,]    8
[1053,]    4
[1054,]    4
[1055,]    1
[1056,]    3
[1057,]    6
[1058,]    4
[1059,]    3
[1060,]    6
[1061,]    3
[1062,]    8
[1063,]    4
[1064,]    1
[1065,]    5
[1066,]    7
[1067,]    5
[1068,]    3
[1069,]    6
[1070,]    4
[1071,]    5
[1072,]    4
[1073,]    6
[1074,]    5
[1075,]    8
[1076,]    9
[1077,]    9
[1078,]    9
[1079,]    7
[1080,]    9
[1081,]    7
[1082,]    7
[1083,]    4
[1084,]    5
[1085,]    4
[1086,]    5
[1087,]    7
[1088,]    4
[1089,]    4
[1090,]    5
[1091,]    4
[1092,]    5
[1093,]    6
[1094,]    7
[1095,]    6
[1096,]    3
[1097,]    6
[1098,]    2
[1099,]    5
[1100,]    9
[1101,]    9
[1102,]    5
[1103,]    5
[1104,]    6
[1105,]    6
[1106,]    6
[1107,]    4
[1108,]    7
[1109,]    6
[1110,]    6
[1111,]    8
[1112,]    7
[1113,]    9
[1114,]    1
[1115,]    6
[1116,]    5
[1117,]    5
[1118,]    4
[1119,]    5
[1120,]    5
[1121,]    5
[1122,]    4
[1123,]    7
[1124,]    5
[1125,]    6
[1126,]    5
[1127,]    3
[1128,]    4
[1129,]    4
[1130,]    8
[1131,]    7
[1132,]    7
[1133,]    4
[1134,]    7
[1135,]    5
[1136,]    7
[1137,]    4
[1138,]    6
[1139,]    8
[1140,]    3
[1141,]    7
[1142,]    7
[1143,]    6
[1144,]    6
[1145,]    3
[1146,]    5
[1147,]    3
[1148,]    5
[1149,]    6
[1150,]    8
[1151,]    7
[1152,]    6
[1153,]    6
[1154,]    9
[1155,]    9
[1156,]    9
[1157,]    9
[1158,]    3
[1159,]    3
[1160,]    3
[1161,]    3
[1162,]    5
[1163,]    4
[1164,]    7
[1165,]    9
[1166,]    7
[1167,]    7
[1168,]    9
[1169,]    5
[1170,]    3
[1171,]    6
[1172,]    6
[1173,]    1
[1174,]    5
[1175,]    7
[1176,]    7
[1177,]    7
[1178,]    5
[1179,]    7
[1180,]    4
[1181,]    5
[1182,]    5
[1183,]    7
[1184,]    5
[1185,]    3
[1186,]    7
[1187,]    5
[1188,]    3
[1189,]    7
[1190,]    7
[1191,]    4
[1192,]    5
[1193,]    5
[1194,]    6
[1195,]    5
[1196,]    3
[1197,]    8
[1198,]    4
[1199,]    4
[1200,]    6
[1201,]    2
[1202,]    2
[1203,]    7
[1204,]    6
[1205,]    9
[1206,]    5
[1207,]    7
[1208,]    7
[1209,]    7
[1210,]    7
[1211,]    5
[1212,]    7
[1213,]    2
[1214,]    8
[1215,]    4
[1216,]    5
[1217,]    6
[1218,]    7
[1219,]    8
[1220,]    3
[1221,]    7
[1222,]    4
[1223,]    6
[1224,]    5
[1225,]    2
[1226,]    4
[1227,]    7
[1228,]    7
[1229,]    8
[1230,]    6
[1231,]    8
[1232,]    8
[1233,]    7
[1234,]    7
[1235,]    5
[1236,]    6
[1237,]    5
[1238,]    3
[1239,]    5
[1240,]    7
[1241,]    6
[1242,]    8
[1243,]    7
[1244,]    7
[1245,]    6
[1246,]    4
[1247,]    5
[1248,]    5
[1249,]    9
[1250,]    2
[1251,]    3
[1252,]    4
[1253,]    4
[1254,]    3
[1255,]    3
[1256,]    5
[1257,]    4
[1258,]    3
[1259,]    5
[1260,]    4
[1261,]    5
[1262,]    8
[1263,]    3
[1264,]    9
[1265,]    5
[1266,]    7
[1267,]    4
[1268,]    5
[1269,]    6
[1270,]    7
[1271,]    7
[1272,]    6
[1273,]    6
[1274,]    8
[1275,]    7
[1276,]    5
[1277,]    8
[1278,]    7
[1279,]    4
[1280,]    4
[1281,]    7
[1282,]    8
[1283,]    2
[1284,]    7
[1285,]    9
[1286,]    8
[1287,]    7
[1288,]    4
[1289,]    6
[1290,]    8
[1291,]    6
[1292,]    2
[1293,]    5
[1294,]    7
[1295,]    8
[1296,]    8
[1297,]    5
[1298,]    7
[1299,]    8
[1300,]    8
[1301,]    3
[1302,]    7
[1303,]    6
[1304,]    3
[1305,]    6
[1306,]    5
[1307,]    5
[1308,]    5
[1309,]    6
[1310,]    3
[1311,]    4
[1312,]    5
[1313,]    2
[1314,]    9
[1315,]    6
[1316,]    2
[1317,]    5
[1318,]    8
[1319,]    6
[1320,]    6
[1321,]    6
[1322,]    5
[1323,]    7
[1324,]    6
[1325,]    7
[1326,]    8
[1327,]    8
[1328,]    9
[1329,]    7
[1330,]    6
[1331,]    4
[1332,]    4
[1333,]    6
[1334,]    7
[1335,]    9
[1336,]    6
[1337,]    5
[1338,]    7
[1339,]    8
[1340,]    6
[1341,]    6
[1342,]    7
[1343,]    5
[1344,]    6
[1345,]    3
[1346,]    9
[1347,]    7
[1348,]    9
[1349,]    5
[1350,]    7
[1351,]    7
[1352,]    7
[1353,]    8
[1354,]    7
[1355,]    9
[1356,]    4
[1357,]    2
[1358,]    5
[1359,]    5
[1360,]    4
[1361,]    4
[1362,]    5
[1363,]    5
[1364,]    3
[1365,]    7
[1366,]    4
[1367,]    6
[1368,]    4
[1369,]    3
[1370,]    8
[1371,]    6
[1372,]    2
[1373,]    9
[1374,]    5
[1375,]    8
[1376,]    5
[1377,]    6
[1378,]    4
[1379,]    3
[1380,]    7
[1381,]    6
[1382,]    8
[1383,]    7
[1384,]    5
[1385,]    4
[1386,]    1
[1387,]    2
[1388,]    3
[1389,]    4
[1390,]    7
[1391,]    6
[1392,]    3
[1393,]    5
[1394,]    4
[1395,]    6
[1396,]    7
[1397,]    8
[1398,]    8
[1399,]    2
[1400,]    8
[1401,]    6
[1402,]    4
[1403,]    6
[1404,]    7
[1405,]    3
[1406,]    5
[1407,]    6
[1408,]    9
[1409,]    6
[1410,]    5
[1411,]    2
[1412,]    6
[1413,]    7
[1414,]    5
[1415,]    6
[1416,]    6
[1417,]    4
[1418,]    7
[1419,]    5
[1420,]    8
[1421,]    6
[1422,]    4
[1423,]    7
[1424,]    7
[1425,]    6
[1426,]    7
[1427,]    7
[1428,]    7
[1429,]    5
[1430,]    8
[1431,]    6
[1432,]    6
[1433,]    5
[1434,]    5
[1435,]    5
[1436,]    7
[1437,]    8
[1438,]    8
[1439,]    5
[1440,]    5
[1441,]    6
[1442,]    6
[1443,]    7
[1444,]    5
[1445,]    5
[1446,]    5
[1447,]    4
[1448,]    4
[1449,]    7
[1450,]    4
[1451,]    5
[1452,]    6
[1453,]    8
[1454,]    6
[1455,]    7
[1456,]    5
[1457,]    7
[1458,]    8
[1459,]    6
[1460,]    9
[1461,]    9
[1462,]    7
[1463,]    5
[1464,]    7
[1465,]    8
[1466,]    6
[1467,]    4
[1468,]    3
[1469,]    6
[1470,]    7
[1471,]    8
[1472,]    8
[1473,]    4
[1474,]    5
[1475,]    6
[1476,]    4
[1477,]    4
[1478,]    3
[1479,]    6
[1480,]    7
[1481,]    4
[1482,]    7
[1483,]    3
[1484,]    4
[1485,]    3
[1486,]    5
[1487,]    3
[1488,]    7
[1489,]    7
[1490,]    5
[1491,]    6
[1492,]    7
[1493,]    9
[1494,]    6
[1495,]    8
[1496,]    2
[1497,]    6
[1498,]    7
[1499,]    6
[1500,]    2
[1501,]    3
[1502,]    4
[1503,]    2
[1504,]    2
[1505,]    5
[1506,]    7
[1507,]    6
[1508,]    3
[1509,]    4
[1510,]    3
[1511,]    5
[1512,]    4
[1513,]    6
[1514,]    9
[1515,]    6
[1516,]    5
[1517,]    3
[1518,]    4
[1519,]    3
[1520,]    6
[1521,]    7
[1522,]    5
[1523,]    7
[1524,]    5
[1525,]    3
[1526,]    6
[1527,]    6
[1528,]    4
[1529,]    5
[1530,]    7
[1531,]    4
[1532,]    3
[1533,]    7
[1534,]    7
[1535,]    4
[1536,]    7
[1537,]    9
[1538,]    7
[1539,]    5
[1540,]    5
[1541,]    3
[1542,]    2
[1543,]    5
[1544,]    4
[1545,]    3
[1546,]    4
[1547,]    6
[1548,]    5
[1549,]    7
[1550,]    6
[1551,]    8
[1552,]    7
[1553,]    5
[1554,]    5
[1555,]    5
[1556,]    8
[1557,]    7
[1558,]    5
[1559,]    6
[1560,]    4
[1561,]    8
[1562,]    6
[1563,]    5
[1564,]    5
[1565,]    5
[1566,]    5
[1567,]    6
[1568,]    3
[1569,]    4
[1570,]    4
[1571,]    6
[1572,]    9
[1573,]    4
[1574,]    6
[1575,]    9
[1576,]    6
[1577,]    3
[1578,]    2
[1579,]    2
[1580,]    6
[1581,]    6
[1582,]    6
[1583,]    2
[1584,]    6
[1585,]    6
[1586,]    3
[1587,]    6
[1588,]    6
[1589,]    8
[1590,]    6
[1591,]    7
[1592,]    7
[1593,]    5
[1594,]    7
[1595,]    8
[1596,]    7
[1597,]    4
[1598,]    5
[1599,]    9
[1600,]    7
[1601,]    8
[1602,]    9
[1603,]    9
[1604,]    7
[1605,]    6
[1606,]    7
[1607,]    8
[1608,]    4
[1609,]    7
[1610,]    4
[1611,]    3
[1612,]    6
[1613,]    5
[1614,]    9
[1615,]    3
[1616,]    6
[1617,]    7
[1618,]    5
[1619,]    8
[1620,]    6
[1621,]    7
[1622,]    4
[1623,]    3
[1624,]    3
[1625,]    7
[1626,]    4
[1627,]    6
[1628,]    4
[1629,]    6
[1630,]    3
[1631,]    6
[1632,]    7
[1633,]    7
[1634,]    5
[1635,]    6
[1636,]    7
[1637,]    8
[1638,]    7
[1639,]    6
[1640,]    7
[1641,]    6
[1642,]    8
[1643,]    8
[1644,]    8
[1645,]    8
[1646,]    3
[1647,]    3
[1648,]    6
[1649,]    5
[1650,]    5
[1651,]    4
[1652,]    8
[1653,]    4
[1654,]    7
[1655,]    6
[1656,]    3
[1657,]    3
[1658,]    5
[1659,]    4
[1660,]    6
[1661,]    6
[1662,]    2
[1663,]    8
[1664,]    9
[1665,]    5
[1666,]    5
[1667,]    7
[1668,]    6
[1669,]    4
[1670,]    7
[1671,]    7
[1672,]    4
[1673,]    6
[1674,]    8
[1675,]    4
[1676,]    6
[1677,]    4
[1678,]    4
[1679,]    6
[1680,]    5
[1681,]    2
[1682,]    1
[1683,]    6
[1684,]    5
[1685,]    5
[1686,]    8
[1687,]    7
[1688,]    5
[1689,]    4
[1690,]    3
[1691,]    4
[1692,]    7
[1693,]    8
[1694,]    8
[1695,]    6
[1696,]    4
[1697,]    7
[1698,]    6
[1699,]    5
[1700,]    5
[1701,]    5
[1702,]    5
[1703,]    3
[1704,]    4
[1705,]    7
[1706,]    2
[1707,]    6
[1708,]    8
[1709,]    8
[1710,]    6
[1711,]    4
[1712,]    7
[1713,]    5
[1714,]    5
[1715,]    5
[1716,]    3
[1717,]    4
[1718,]    6
[1719,]    9
[1720,]    7
[1721,]    4
[1722,]    4
[1723,]    7
[1724,]    5
[1725,]    8
[1726,]    3
[1727,]    9
[1728,]    4
[1729,]    5
[1730,]    4
[1731,]    6
[1732,]    8
[1733,]    6
[1734,]    7
[1735,]    4
[1736,]    5
[1737,]    6
[1738,]    5
[1739,]    5
[1740,]    4
[1741,]    5
[1742,]    7
[1743,]    8
[1744,]    2
[1745,]    4
[1746,]    5
[1747,]    5
[1748,]    8
[1749,]    5
[1750,]    7
[1751,]    6
[1752,]    4
[1753,]    7
[1754,]    5
[1755,]    3
[1756,]    6
[1757,]    1
[1758,]    2
[1759,]    2
[1760,]    7
[1761,]    7
[1762,]    8
[1763,]    5
[1764,]    5
[1765,]    5
[1766,]    5
[1767,]    7
[1768,]    8
[1769,]    5
[1770,]    7
[1771,]    8
[1772,]    7
[1773,]    7
[1774,]    7
[1775,]    5
[1776,]    7
[1777,]    6
[1778,]    5
[1779,]    4
[1780,]    4
[1781,]    3
[1782,]    5
[1783,]    7
[1784,]    8
[1785,]    6
[1786,]    5
[1787,]    6
[1788,]    8
[1789,]    3
[1790,]    5
[1791,]    4
[1792,]    5
[1793,]    5
[1794,]    5
[1795,]    7
[1796,]    6
[1797,]    4
[1798,]    4
[1799,]    6
[1800,]    3
[1801,]    2
[1802,]    5
[1803,]    7
[1804,]    8
[1805,]    6
[1806,]    8
[1807,]    7
[1808,]    9
[1809,]    8
[1810,]    7
[1811,]    9
[1812,]    8
[1813,]    3
[1814,]    8
[1815,]    7
[1816,]    6
[1817,]    2
[1818,]    5
[1819,]    1
[1820,]    7
[1821,]    6
[1822,]    4
[1823,]    8
[1824,]    8
[1825,]    8
[1826,]    4
[1827,]    7
[1828,]    6
[1829,]    7
[1830,]    6
[1831,]    7
[1832,]    7
[1833,]    6
[1834,]    9
[1835,]    7
[1836,]    5
[1837,]    8
[1838,]    6
[1839,]    7
[1840,]    7
[1841,]    7
[1842,]    6
[1843,]    8
[1844,]    4
[1845,]    8
[1846,]    5
[1847,]    4
[1848,]    3
[1849,]    7
[1850,]    8
[1851,]    7
[1852,]    8
[1853,]    5
[1854,]    4
[1855,]    5
[1856,]    7
[1857,]    7
[1858,]    7
[1859,]    8
[1860,]    6
[1861,]    7
[1862,]    5
[1863,]    8
[1864,]    5
[1865,]    5
[1866,]    7
[1867,]    8
[1868,]    4
[1869,]    8
[1870,]    8
[1871,]    7
[1872,]    9
[1873,]    5
[1874,]    9
[1875,]    9
[1876,]    7
[1877,]    5
[1878,]    7
[1879,]    5
[1880,]    5
[1881,]    5
[1882,]    4
[1883,]    4
[1884,]    5
[1885,]    8
[1886,]    4
[1887,]    5
[1888,]    3
[1889,]    4
[1890,]    5
[1891,]    7
[1892,]    8
[1893,]    4
[1894,]    9
[1895,]    5
[1896,]    7
[1897,]    3
[1898,]    5
[1899,]    4
[1900,]    5
[1901,]    6
[1902,]    7
[1903,]    2
[1904,]    9
[1905,]    7
[1906,]    7
[1907,]    6
[1908,]    9
[1909,]    6
[1910,]    8
[1911,]    7
[1912,]    3
[1913,]    7
[1914,]    8
[1915,]    9
[1916,]    6
[1917,]    7
[1918,]    3
[1919,]    5
[1920,]    8
[1921,]    5
[1922,]    6
[1923,]    7
[1924,]    8
[1925,]    7
[1926,]    6
[1927,]    5
[1928,]    8
[1929,]    5
[1930,]    6
[1931,]    8
[1932,]    8
[1933,]    3
[1934,]    6
[1935,]    8
[1936,]    5
[1937,]    9
[1938,]    5
[1939,]    0
[1940,]    3
[1941,]    5
[1942,]    7
[1943,]    4
[1944,]    3
[1945,]    5
[1946,]    4
[1947,]    5
[1948,]    7
[1949,]    5
[1950,]    9
[1951,]    9
[1952,]    8
[1953,]    9
[1954,]    6
[1955,]    3
[1956,]    5
[1957,]    9
[1958,]    8
[1959,]    8
[1960,]    9
[1961,]    8
[1962,]    8
[1963,]    7
[1964,]    6
[1965,]    5
[1966,]    6
[1967,]    4
[1968,]    6
[1969,]    8
[1970,]    8
[1971,]    6
[1972,]    5
[1973,]    8
[1974,]    5
[1975,]    7
[1976,]    7
[1977,]    6
[1978,]    1
[1979,]    3
[1980,]    6
[1981,]    6
[1982,]    7
[1983,]    5
[1984,]    6
[1985,]    6
[1986,]    4
[1987,]    8
[1988,]    3
[1989,]    5
[1990,]    6
[1991,]    5
[1992,]    4
[1993,]    6
[1994,]    6
[1995,]    8
[1996,]    7
[1997,]    5
[1998,]    9
[1999,]    5
[2000,]    9
[2001,]    5
[2002,]    6
[2003,]    2
[2004,]    3
[2005,]    8
[2006,]    6
[2007,]    6
[2008,]    7
[2009,]    6
[2010,]    8
[2011,]    9
[2012,]    8
[2013,]    8
[2014,]    9
[2015,]    6
[2016,]    6
[2017,]    3
[2018,]    6
[2019,]    4
[2020,]    5
[2021,]    6
[2022,]    6
[2023,]    8
[2024,]    6
[2025,]    4
[2026,]    7
[2027,]    7
[2028,]    9
[2029,]    6
[2030,]    7
[2031,]    5
[2032,]    6
[2033,]    6
[2034,]    7
[2035,]    7
[2036,]    7
[2037,]    7
[2038,]    8
[2039,]    5
[2040,]    3
[2041,]    6
[2042,]    6
[2043,]    4
[2044,]    5
[2045,]    7
[2046,]    5
[2047,]    5
[2048,]    4
[2049,]    6
[2050,]    8
[2051,]    4
[2052,]    8
[2053,]    7
[2054,]    4
[2055,]    4
[2056,]    7
[2057,]    7
[2058,]    4
[2059,]    4
[2060,]    8
[2061,]    7
[2062,]    7
[2063,]    5
[2064,]    8
[2065,]    3
[2066,]    5
[2067,]    5
[2068,]    5
[2069,]    8
[2070,]    6
[2071,]    6
[2072,]    7
[2073,]    6
[2074,]    2
[2075,]    3
[2076,]    8
[2077,]    8
[2078,]    5
[2079,]    4
[2080,]    5
[2081,]    3
[2082,]    5
[2083,]    4
[2084,]    5
[2085,]    4
[2086,]    7
[2087,]    8
[2088,]    7
[2089,]    7
[2090,]    5
[2091,]    7
[2092,]    2
[2093,]    1
[2094,]    1
[2095,]    4
[2096,]    4
[2097,]    6
[2098,]    7
[2099,]    6
[2100,]    8
[2101,]    8
[2102,]    8
[2103,]    7
[2104,]    8
[2105,]    6
[2106,]    4
[2107,]    1
[2108,]    5
[2109,]    6
[2110,]    7
[2111,]    5
[2112,]    8
[2113,]    4
[2114,]    5
[2115,]    8
[2116,]    6
[2117,]    5
[2118,]    3
[2119,]    4
[2120,]    5
[2121,]    7
[2122,]    8
[2123,]    6
[2124,]    6
[2125,]    9
[2126,]    7
[2127,]    8
[2128,]    7
[2129,]    6
[2130,]    5
[2131,]    6
[2132,]    5
[2133,]    5
[2134,]    5
[2135,]    7
[2136,]    5
[2137,]    7
[2138,]    3
[2139,]    6
[2140,]    4
[2141,]    4
[2142,]    5
[2143,]    5
[2144,]    5
[2145,]    7
[2146,]    9
[2147,]    5
[2148,]    3
[2149,]    3
[2150,]    5
[2151,]    4
[2152,]    5
[2153,]    2
[2154,]    9
[2155,]    7
[2156,]    5
[2157,]    6
[2158,]    3
[2159,]    4
[2160,]    6
[2161,]    4
[2162,]    7
[2163,]    3
[2164,]    7
[2165,]    6
[2166,]    8
[2167,]    2
[2168,]    9
[2169,]    5
[2170,]    5
[2171,]    6
[2172,]    4
[2173,]    2
[2174,]    7
[2175,]    4
[2176,]    6
[2177,]    9
[2178,]    6
[2179,]    8
[2180,]    7
[2181,]    9
[2182,]    4
[2183,]    2
[2184,]    1
[2185,]    3
[2186,]    4
[2187,]    3
[2188,]    5
[2189,]    5
[2190,]    2
[2191,]    2
[2192,]    5
[2193,]    2
[2194,]    3
[2195,]    6
[2196,]    4
[2197,]    6
[2198,]    7
[2199,]    5
[2200,]    4
[2201,]    4
[2202,]    3
[2203,]    8
[2204,]    3
[2205,]    6
[2206,]    7
[2207,]    5
[2208,]    6
[2209,]    7
[2210,]    5
[2211,]    5
[2212,]    7
[2213,]    4
[2214,]    8
[2215,]    4
[2216,]    8
[2217,]    4
[2218,]    5
[2219,]    8
[2220,]    3
[2221,]    5
[2222,]    8
[2223,]    6
[2224,]    5
[2225,]    5
[2226,]    9
[2227,]    8
[2228,]    9
[2229,]    7
[2230,]    8
[2231,]    4
[2232,]    6
[2233,]    5
[2234,]    4
[2235,]    5
[2236,]    7
[2237,]    4
[2238,]    8
[2239,]    4
[2240,]    4
[2241,]    6
[2242,]    7
[2243,]    6
[2244,]    8
[2245,]    7
[2246,]    5
[2247,]    8
[2248,]    6
[2249,]    6
[2250,]    8
[2251,]    6
[2252,]    3
[2253,]    2
[2254,]    7
[2255,]    8
[2256,]    9
[2257,]    8
[2258,]    2
[2259,]    7
[2260,]    7
[2261,]    8
[2262,]    5
[2263,]    5
[2264,]    5
[2265,]    7
[2266,]    6
[2267,]    6
[2268,]    7
[2269,]    5
[2270,]    5
[2271,]    6
[2272,]    4
[2273,]    7
[2274,]    5
[2275,]    3
[2276,]    7
[2277,]    6
[2278,]    5
[2279,]    6
[2280,]    5
[2281,]    6
[2282,]    7
[2283,]    6
[2284,]    6
[2285,]    7
[2286,]    5
[2287,]    6
[2288,]    5
[2289,]    4
[2290,]    6
[2291,]    5
[2292,]    5
[2293,]    6
[2294,]    7
[2295,]    5
[2296,]    4
[2297,]    7
[2298,]    4
[2299,]    3
[2300,]    2
[2301,]    6
[2302,]    7
[2303,]    4
[2304,]    2
[2305,]    6
[2306,]    7
[2307,]    4
[2308,]    8
[2309,]    5
[2310,]    6
[2311,]    7
[2312,]    5
[2313,]    4
[2314,]    7
[2315,]    6
[2316,]    7
[2317,]    8
[2318,]    6
[2319,]    5
[2320,]    0
[2321,]    9
[2322,]    6
[2323,]    5
[2324,]    6
[2325,]    7
[2326,]    8
[2327,]    3
[2328,]    8
[2329,]    6
[2330,]    4
[2331,]    6
[2332,]    8
[2333,]    3
[2334,]    7
[2335,]    3
[2336,]    5
[2337,]    7
[2338,]    7
[2339,]    8
[2340,]    8
[2341,]    9
[2342,]    6
[2343,]    3
[2344,]    4
[2345,]    9
[2346,]    8
[2347,]    7
[2348,]    8
[2349,]    6
[2350,]    5
[2351,]    6
[2352,]    7
[2353,]    6
[2354,]    5
[2355,]    1
[2356,]    2
[2357,]    4
[2358,]    2
[2359,]    6
[2360,]    8
[2361,]    7
[2362,]    3
[2363,]    9
[2364,]    5
[2365,]    6
[2366,]    4
[2367,]    7
[2368,]    4
[2369,]    8
[2370,]    8
[2371,]    6
[2372,]    4
[2373,]    4
[2374,]    4
[2375,]    9
[2376,]    8
[2377,]    4
[2378,]    7
[2379,]    5
[2380,]    7
[2381,]    6
[2382,]    5
[2383,]    7
[2384,]    6
[2385,]    5
[2386,]    6
[2387,]    5
[2388,]    6
[2389,]    8
[2390,]    5
[2391,]    8
[2392,]    8
[2393,]    6
[2394,]    2
[2395,]    5
[2396,]    5
[2397,]    6
[2398,]    3
[2399,]    5
[2400,]    8
[2401,]    8
[2402,]    7
[2403,]    6
[2404,]    8
[2405,]    6
[2406,]    5
[2407,]    7
[2408,]    7
[2409,]    7
[2410,]    8
[2411,]    8
[2412,]    7
[2413,]    2
[2414,]    5
[2415,]    6
[2416,]    5
[2417,]    7
[2418,]    6
[2419,]    6
[2420,]    5
[2421,]    6
[2422,]    6
[2423,]    7
[2424,]    3
[2425,]    5
[2426,]    3
[2427,]    6
[2428,]    7
[2429,]    7
[2430,]    7
[2431,]    3
[2432,]    7
[2433,]    8
[2434,]    7
[2435,]    6
[2436,]    5
[2437,]    4
[2438,]    6
[2439,]    8
[2440,]    4
[2441,]    7
[2442,]    4
[2443,]    5
[2444,]    4
[2445,]    6
[2446,]    8
[2447,]    5
[2448,]    5
[2449,]    5
[2450,]    6
[2451,]    5
[2452,]    9
[2453,]    7
[2454,]    3
[2455,]    8
[2456,]    6
[2457,]    6
[2458,]    7
[2459,]    6
[2460,]    5
[2461,]    6
[2462,]    7
[2463,]    2
[2464,]    7
[2465,]    8
[2466,]    6
[2467,]    8
[2468,]    8
[2469,]    3
[2470,]    7
[2471,]    8
[2472,]    6
[2473,]    4
[2474,]    7
[2475,]    5
[2476,]    6
[2477,]    8
[2478,]    7
[2479,]    6
[2480,]    7
[2481,]    5
[2482,]    4
[2483,]    6
[2484,]    9
[2485,]    7
[2486,]    7
[2487,]    9
[2488,]    6
[2489,]    7
[2490,]    7
[2491,]    8
[2492,]    9
[2493,]    5
[2494,]    5
[2495,]    2
[2496,]    4
[2497,]    5
[2498,]    4
[2499,]    4
[2500,]    1
[2501,]    6
[2502,]    6
[2503,]    5
[2504,]    5
[2505,]    9
[2506,]    4
[2507,]    5
[2508,]    4
[2509,]    3
[2510,]    8
[2511,]    3
[2512,]    5
[2513,]    6
[2514,]    7
[2515,]    6
[2516,]    4
[2517,]    9
[2518,]    6
[2519,]    6
[2520,]    5
[2521,]    6
[2522,]    5
[2523,]    6
[2524,]    5
[2525,]    4
[2526,]    6
[2527,]    8
[2528,]    4
[2529,]    4
[2530,]    2
[2531,]    7
[2532,]    7
[2533,]    9
[2534,]    6
[2535,]    9
[2536,]    5
[2537,]    7
[2538,]    8
[2539,]    3
[2540,]    4
[2541,]    5
[2542,]    5
[2543,]    8
[2544,]    7
[2545,]    8
[2546,]    4
[2547,]    4
[2548,]    4
[2549,]    6
[2550,]    8
[2551,]    6
[2552,]    4
[2553,]    6
[2554,]    8
[2555,]    3
[2556,]    8
[2557,]    6
[2558,]    6
[2559,]    4
[2560,]    6
[2561,]    5
[2562,]    3
[2563,]    6
[2564,]    9
[2565,]    5
[2566,]    7
[2567,]    3
[2568,]    3
[2569,]    4
[2570,]    8
[2571,]    7
[2572,]    8
[2573,]    4
[2574,]    6
[2575,]    7
[2576,]    6
[2577,]    3
[2578,]    5
[2579,]    9
[2580,]    7
[2581,]    6
[2582,]    5
[2583,]    6
[2584,]    6
[2585,]    7
[2586,]    8
[2587,]    9
[2588,]    6
[2589,]    9
[2590,]    9
[2591,]    4
[2592,]    6
[2593,]    5
[2594,]    3
[2595,]    4
[2596,]    4
[2597,]    5
[2598,]    8
[2599,]    2
[2600,]    9
[2601,]    6
[2602,]    8
[2603,]    7
[2604,]    7
[2605,]    6
[2606,]    7
[2607,]    4
[2608,]    4
[2609,]    5
[2610,]    8
[2611,]    5
[2612,]    8
[2613,]    8
[2614,]    7
[2615,]    6
[2616,]    6
[2617,]    2
[2618,]    3
[2619,]    5
[2620,]    6
[2621,]    5
[2622,]    5
[2623,]    6
[2624,]    7
[2625,]    6
[2626,]    5
[2627,]    5
[2628,]    7
[2629,]    5
[2630,]    7
[2631,]    8
[2632,]    8
[2633,]    8
[2634,]    8
[2635,]    8
[2636,]    7
[2637,]    8
[2638,]    9
[2639,]    8
[2640,]    7
[2641,]    7
[2642,]    7
[2643,]    6
[2644,]    8
[2645,]    7
[2646,]    4
[2647,]    6
[2648,]    5
[2649,]    5
[2650,]    7
[2651,]    3
[2652,]    5
[2653,]    3
[2654,]    5
[2655,]    7
[2656,]    5
[2657,]    7
[2658,]    7
[2659,]    6
[2660,]    7
[2661,]    6
[2662,]    8
[2663,]    6
[2664,]    5
[2665,]    6
[2666,]    5
[2667,]    7
[2668,]    6
[2669,]    7
[2670,]    5
[2671,]    4
[2672,]    5
[2673,]    7
[2674,]    6
[2675,]    7
[2676,]    7
[2677,]    4
[2678,]    9
[2679,]    3
[2680,]    4
[2681,]    7
[2682,]    6
[2683,]    7
[2684,]    6
[2685,]    5
[2686,]    6
[2687,]    7
[2688,]    5
[2689,]    9
[2690,]    6
[2691,]    4
[2692,]    8
[2693,]    2
[2694,]    6
[2695,]    5
[2696,]    7
[2697,]    7
[2698,]    6
[2699,]    9
[2700,]    6
[2701,]    6
[2702,]    7
[2703,]    6
[2704,]    6
[2705,]    7
[2706,]    6
[2707,]    3
[2708,]    3
[2709,]    7
[2710,]    9
[2711,]    7
[2712,]    8
[2713,]    8
[2714,]    4
[2715,]    6
[2716,]    7
[2717,]    9
[2718,]    8
[2719,]    1
[2720,]    4
[2721,]    5
[2722,]    6
[2723,]    3
[2724,]    6
[2725,]    6
[2726,]    7
[2727,]    5
[2728,]    6
[2729,]    3
[2730,]    5
[2731,]    3
[2732,]    7
[2733,]    7
[2734,]    7
[2735,]    6
[2736,]    7
[2737,]    8
[2738,]    9
[2739,]    4
[2740,]    6
[2741,]    7
[2742,]    7
[2743,]    5
[2744,]    5
[2745,]    7
[2746,]    8
[2747,]    7
[2748,]    6
[2749,]    3
[2750,]    5
[2751,]    7
[2752,]    4
[2753,]    6
[2754,]    5
[2755,]    7
[2756,]    7
[2757,]    4
[2758,]    7
[2759,]    6
[2760,]    6
[2761,]    5
[2762,]    4
[2763,]    6
[2764,]    7
[2765,]    6
[2766,]    5
[2767,]    5
[2768,]    3
[2769,]    7
[2770,]    8
[2771,]    9
[2772,]    8
[2773,]    8
[2774,]    7
[2775,]    7
[2776,]    4
[2777,]    5
[2778,]    3
[2779,]    6
[2780,]    6
[2781,]    7
[2782,]    7
[2783,]    6
[2784,]    5
[2785,]    4
[2786,]    6
[2787,]    6
[2788,]    7
[2789,]    8
[2790,]    7
[2791,]    4
[2792,]    5
[2793,]    8
[2794,]    9
[2795,]    1
[2796,]    7
[2797,]    6
[2798,]    7
[2799,]    1
[2800,]    7
[2801,]    7
[2802,]    8
[2803,]    8
[2804,]    7
[2805,]    8
[2806,]    7
[2807,]    7
[2808,]    7
[2809,]    5
[2810,]    8
[2811,]    6
[2812,]    5
[2813,]    8
[2814,]    8
[2815,]    7
[2816,]    5
[2817,]    7
[2818,]    5
[2819,]    3
[2820,]    7
[2821,]    7
[2822,]    4
[2823,]    6
[2824,]    2
[2825,]    6
[2826,]    5
[2827,]    5
[2828,]    3
[2829,]    3
[2830,]    3
[2831,]    5
[2832,]    5
[2833,]    6
[2834,]    5
[2835,]    6
[2836,]    7
[2837,]    9
[2838,]    7
[2839,]    5
[2840,]    7
[2841,]    8
[2842,]    5
[2843,]    8
[2844,]    5
[2845,]    5
[2846,]    8
[2847,]    7
[2848,]    7
[2849,]    0
[2850,]    4
[2851,]    7
[2852,]    6
[2853,]    6
[2854,]    7
[2855,]    7
[2856,]    6
[2857,]    3
[2858,]    2
[2859,]    6
[2860,]    6
[2861,]    2
[2862,]    3
[2863,]    3
[2864,]    2
[2865,]    6
[2866,]    7
[2867,]    2
[2868,]    4
[2869,]    6
[2870,]    6
[2871,]    9
[2872,]    5
[2873,]    6
[2874,]    7
[2875,]    4
[2876,]    5
[2877,]    3
[2878,]    5
[2879,]    6
[2880,]    5
[2881,]    8
[2882,]    8
[2883,]    4
[2884,]    4
[2885,]    5
[2886,]    6
[2887,]    5
[2888,]    9
[2889,]    7
[2890,]    8
[2891,]    9
[2892,]    6
[2893,]    7
[2894,]    7
[2895,]    7
[2896,]    5
[2897,]    5
[2898,]    8
[2899,]    5
[2900,]    5
[2901,]    6
[2902,]    4
[2903,]    7
[2904,]    6
[2905,]    6
[2906,]    5
[2907,]    8
[2908,]    6
[2909,]    5
[2910,]    7
[2911,]    6
[2912,]    5
[2913,]    5
[2914,]    6
[2915,]    5
[2916,]    7
[2917,]    7
[2918,]    4
[2919,]    4
[2920,]    6
[2921,]    5
[2922,]    7
[2923,]    5
[2924,]    4
[2925,]    7
[2926,]    3
[2927,]    6
[2928,]    7
[2929,]    2
[2930,]    4
[2931,]    6
[2932,]    6
[2933,]    5
[2934,]    5
[2935,]    5
[2936,]    8
[2937,]    5
[2938,]    3
[2939,]    6
[2940,]    8
[2941,]    7
[2942,]    7
[2943,]    8
[2944,]    6
[2945,]    6
[2946,]    4
[2947,]    5
[2948,]    9
[2949,]    7
[2950,]    4
[2951,]    4
[2952,]    5
[2953,]    4
[2954,]    9
[2955,]    5
[2956,]    5
[2957,]    1
[2958,]    3
[2959,]    5
[2960,]    7
[2961,]    7
[2962,]    4
[2963,]    8
[2964,]    6
[2965,]    4
[2966,]    3
[2967,]    6
[2968,]    7
[2969,]    3
[2970,]    3
[2971,]    4
[2972,]    5
[2973,]    6
[2974,]    7
[2975,]    4
[2976,]    7
[2977,]    5
[2978,]    7
[2979,]    5
[2980,]    7
[2981,]    6
[2982,]    7
[2983,]    6
[2984,]    7
[2985,]    5
[2986,]    6
[2987,]    6
[2988,]    7
[2989,]    7
[2990,]    5
[2991,]    8
[2992,]    3
[2993,]    4
[2994,]    4
[2995,]    2
[2996,]    5
[2997,]    7
[2998,]    7
[2999,]    6
[3000,]    7
[3001,]    6
[3002,]    5
[3003,]    7
[3004,]    5
[3005,]    6
[3006,]    8
[3007,]    8
[3008,]    7
[3009,]    4
[3010,]    4
[3011,]    4
[3012,]    8
[3013,]    6
[3014,]    5
[3015,]    3
[3016,]    7
[3017,]    8
[3018,]    5
[3019,]    7
[3020,]    6
[3021,]    8
[3022,]    7
[3023,]    6
[3024,]    5
[3025,]    4
[3026,]    4
[3027,]    4
[3028,]    3
[3029,]    7
[3030,]    6
[3031,]    8
[3032,]    6
[3033,]    5
[3034,]    4
[3035,]    4
[3036,]    6
[3037,]    7
[3038,]    8
[3039,]    3
[3040,]    5
[3041,]    5
[3042,]    6
[3043,]    7
[3044,]    6
[3045,]    4
[3046,]    8
[3047,]    6
[3048,]    6
[3049,]    7
[3050,]    8
[3051,]    6
[3052,]    7
[3053,]    7
[3054,]    6
[3055,]    6
[3056,]    5
[3057,]    3
[3058,]    7
[3059,]    7
[3060,]    5
[3061,]    5
[3062,]    7
[3063,]    3
[3064,]    7
[3065,]    7
[3066,]    8
[3067,]    6
[3068,]    4
[3069,]    5
[3070,]    2
[3071,]    7
[3072,]    8
[3073,]    6
[3074,]    7
[3075,]    6
[3076,]    6
[3077,]    7
[3078,]    5
[3079,]    4
[3080,]    3
[3081,]    0
[3082,]    5
[3083,]    6
[3084,]    7
[3085,]    5
[3086,]    5
[3087,]    4
[3088,]    5
[3089,]    5
[3090,]    7
[3091,]    5
[3092,]    2
[3093,]    3
[3094,]    4
[3095,]    7
[3096,]    4
[3097,]    6
[3098,]    4
[3099,]    6
[3100,]    7
[3101,]    4
[3102,]    7
[3103,]    8
[3104,]    8
[3105,]    5
[3106,]    7
[3107,]    8
[3108,]    7
[3109,]    9
[3110,]    7
[3111,]    7
[3112,]    6
[3113,]    7
[3114,]    2
[3115,]    5
[3116,]    7
[3117,]    8
[3118,]    6
[3119,]    6
[3120,]    3
[3121,]    3
[3122,]    5
[3123,]    5
[3124,]    5
[3125,]    6
[3126,]    6
[3127,]    3
[3128,]    4
[3129,]    8
[3130,]    9
[3131,]    9
[3132,]    6
[3133,]    4
[3134,]    3
[3135,]    3
[3136,]    7
[3137,]    5
[3138,]    3
[3139,]    5
[3140,]    6
[3141,]    7
[3142,]    7
[3143,]    6
[3144,]    9
[3145,]    7
[3146,]    7
[3147,]    8
[3148,]    6
[3149,]    7
[3150,]    5
[3151,]    9
[3152,]    7
[3153,]    7
[3154,]    6
[3155,]    6
[3156,]    3
[3157,]    6
[3158,]    7
[3159,]    6
[3160,]    7
[3161,]    6
[3162,]    7
[3163,]    9
[3164,]    6
[3165,]    5
[3166,]    7
[3167,]    5
[3168,]    7
[3169,]    7
[3170,]    6
[3171,]    7
[3172,]    2
[3173,]    9
[3174,]    6
[3175,]    5
[3176,]    8
[3177,]    5
[3178,]    8
[3179,]    5
[3180,]    7
[3181,]    4
[3182,]    6
[3183,]    6
[3184,]    7
[3185,]    3
[3186,]    8
[3187,]    5
[3188,]    8
[3189,]    9
[3190,]    7
[3191,]    8
[3192,]    4
[3193,]    6
[3194,]    2
[3195,]    9
[3196,]    6
[3197,]    8
[3198,]    4
[3199,]    4
[3200,]    7
[3201,]    4
[3202,]    1
[3203,]    6
[3204,]    6
[3205,]    4
[3206,]    6
[3207,]    5
[3208,]    6
[3209,]    5
[3210,]    8
[3211,]    5
[3212,]    3
[3213,]    6
[3214,]    7
[3215,]    6
[3216,]    4
[3217,]    1
[3218,]    3
[3219,]    3
[3220,]    5
[3221,]    2
[3222,]    5
[3223,]    8
[3224,]    6
[3225,]    5
[3226,]    8
[3227,]    5
[3228,]    7
[3229,]    5
[3230,]    7
[3231,]    2
[3232,]    4
[3233,]    4
[3234,]    5
[3235,]    7
[3236,]    5
[3237,]    6
[3238,]    7
[3239,]    8
[3240,]    9
[3241,]    6
[3242,]    4
[3243,]    9
[3244,]    6
[3245,]    5
[3246,]    5
[3247,]    3
[3248,]    1
[3249,]    3
[3250,]    3
[3251,]    5
[3252,]    6
[3253,]    8
[3254,]    6
[3255,]    8
[3256,]    6
[3257,]    6
[3258,]    8
[3259,]    7
[3260,]    7
[3261,]    6
[3262,]    8
[3263,]    4
[3264,]    7
[3265,]    7
[3266,]    6
[3267,]    8
[3268,]    9
[3269,]    8
[3270,]    5
[3271,]    6
[3272,]    6
[3273,]    5
[3274,]    3
[3275,]    1
[3276,]    4
[3277,]    3
[3278,]    8
[3279,]    5
[3280,]    5
[3281,]    6
[3282,]    3
[3283,]    7
[3284,]    8
[3285,]    5
[3286,]    8
[3287,]    5
[3288,]    6
[3289,]    4
[3290,]    4
[3291,]    3
[3292,]    3
[3293,]    3
[3294,]    6
[3295,]    7
[3296,]    7
[3297,]    9
[3298,]    6
[3299,]    5
[3300,]    6
[3301,]    5
[3302,]    5
[3303,]    5
[3304,]    9
[3305,]    7
[3306,]    6
[3307,]    6
[3308,]    7
[3309,]    8
[3310,]    6
[3311,]    6
[3312,]    4
[3313,]    8
[3314,]    5
[3315,]    3
[3316,]    3
[3317,]    2
[3318,]    7
[3319,]    3
[3320,]    3
[3321,]    4
[3322,]    4
[3323,]    2
[3324,]    4
[3325,]    7
[3326,]    6
[3327,]    3
[3328,]    6
[3329,]    3
[3330,]    4
[3331,]    8
[3332,]    5
[3333,]    7
[3334,]    6
[3335,]    5
[3336,]    8
[3337,]    8
[3338,]    7
[3339,]    6
[3340,]    7
[3341,]    5
[3342,]    7
[3343,]    7
[3344,]    9
[3345,]    2
[3346,]    6
[3347,]    5
[3348,]    8
[3349,]    9
[3350,]    6
[3351,]    7
[3352,]    9
[3353,]    3
[3354,]    6
[3355,]    5
[3356,]    7
[3357,]    6
[3358,]    3
[3359,]    5
[3360,]    3
[3361,]    9
[3362,]    7
[3363,]    3
[3364,]    2
[3365,]    5
[3366,]    7
[3367,]    7
[3368,]    8
[3369,]    4
[3370,]    3
[3371,]    0
[3372,]    4
[3373,]    5
[3374,]    7
[3375,]    8
[3376,]    6
[3377,]    4
[3378,]    3
[3379,]    9
[3380,]    9
[3381,]    5
[3382,]    6
[3383,]    3
[3384,]    7
[3385,]    5
[3386,]    5
[3387,]    7
[3388,]    6
[3389,]    7
[3390,]    9
[3391,]    7
[3392,]    8
[3393,]    7
[3394,]    9
[3395,]    6
[3396,]    6
[3397,]    7
[3398,]    5
[3399,]    8
[3400,]    6
[3401,]    2
[3402,]    9
[3403,]    4
[3404,]    3
[3405,]    4
[3406,]    8
[3407,]    8
[3408,]    7
[3409,]    8
[3410,]    9
[3411,]    7
[3412,]    7
[3413,]    4
[3414,]    4
[3415,]    9
[3416,]    7
[3417,]    7
[3418,]    6
[3419,]    9
[3420,]    7
[3421,]    8
[3422,]    3
[3423,]    2
[3424,]    7
[3425,]    7
[3426,]    8
[3427,]    7
[3428,]    7
[3429,]    8
[3430,]    0
[3431,]    7
[3432,]    6
[3433,]    9
[3434,]    6
[3435,]    5
[3436,]    4
[3437,]    7
[3438,]    6
[3439,]    6
[3440,]    6
[3441,]    6
[3442,]    8
[3443,]    6
[3444,]    7
[3445,]    3
[3446,]    5
[3447,]    6
[3448,]    4
[3449,]    8
[3450,]    7
[3451,]    5
[3452,]    9
[3453,]    8
[3454,]    7
[3455,]    7
[3456,]    5
[3457,]    8
[3458,]    7
[3459,]    8
[3460,]    8
[3461,]    6
[3462,]    8
[3463,]    6
[3464,]    8
[3465,]    4
[3466,]    6
[3467,]    7
[3468,]    8
[3469,]    7
[3470,]    8
[3471,]    7
[3472,]    5
[3473,]    8
[3474,]    4
[3475,]    8
[3476,]    5
[3477,]    8
[3478,]    6
[3479,]    5
[3480,]    5
[3481,]    4
[3482,]    5
[3483,]    5
[3484,]    5
[3485,]    8
[3486,]    8
[3487,]    6
[3488,]    4
[3489,]    2
[3490,]    7
[3491,]    8
[3492,]    4
[3493,]    3
[3494,]    5
[3495,]    8
[3496,]    8
[3497,]    8
[3498,]    8
[3499,]    5
[3500,]    2
[3501,]    6
[3502,]    4
[3503,]    8
[3504,]    8
[3505,]    6
[3506,]    6
[3507,]    7
[3508,]    4
[3509,]    3
[3510,]    3
[3511,]    3
[3512,]    7
[3513,]    5
[3514,]    6
[3515,]    9
[3516,]    8
[3517,]    1
[3518,]    7
[3519,]    7
[3520,]    5
[3521,]    5
[3522,]    6
[3523,]    6
[3524,]    3
[3525,]    5
[3526,]    7
[3527,]    9
[3528,]    7
[3529,]    4
[3530,]    4
[3531,]    5
[3532,]    6
[3533,]    5
[3534,]    8
[3535,]    5
[3536,]    7
[3537,]    5
[3538,]    6
[3539,]    4
[3540,]    5
[3541,]    5
[3542,]    7
[3543,]    5
[3544,]    9
[3545,]    5
[3546,]    7
[3547,]    5
[3548,]    8
[3549,]    8
[3550,]    8
[3551,]    3
[3552,]    5
[3553,]    8
[3554,]    5
[3555,]    6
[3556,]    7
[3557,]    4
[3558,]    5
[3559,]    5
[3560,]    2
[3561,]    8
[3562,]    7
[3563,]    7
[3564,]    6
[3565,]    4
[3566,]    7
[3567,]    5
[3568,]    3
[3569,]    5
[3570,]    6
[3571,]    1
[3572,]    8
[3573,]    7
[3574,]    9
[3575,]    5
[3576,]    6
[3577,]    8
[3578,]    5
[3579,]    7
[3580,]    5
[3581,]    5
[3582,]    7
[3583,]    7
[3584,]    8
[3585,]    6
[3586,]    6
[3587,]    5
[3588,]    3
[3589,]    5
[3590,]    7
[3591,]    6
[3592,]    5
[3593,]    6
[3594,]    7
[3595,]    6
[3596,]    7
[3597,]    8
[3598,]    4
[3599,]    5
[3600,]    7
[3601,]    6
[3602,]    5
[3603,]    6
[3604,]    6
[3605,]    8
[3606,]    4
[3607,]    6
[3608,]    4
[3609,]    3
[3610,]    7
[3611,]    3
[3612,]    6
[3613,]    5
[3614,]    6
[3615,]    6
[3616,]    5
[3617,]    6
[3618,]    6
[3619,]    7
[3620,]    5
[3621,]    6
[3622,]    9
[3623,]    6
[3624,]    6
[3625,]    7
[3626,]    4
[3627,]    4
[3628,]    4
[3629,]    2
[3630,]    3
[3631,]    7
[3632,]    9
[3633,]    3
[3634,]    4
[3635,]    8
[3636,]    3
[3637,]    3
[3638,]    7
[3639,]    7
[3640,]    8
[3641,]    5
[3642,]    6
[3643,]    5
[3644,]    6
[3645,]    5
[3646,]    5
[3647,]    4
[3648,]    5
[3649,]    9
[3650,]    9
[3651,]    9
[3652,]    8
[3653,]    8
[3654,]    9
[3655,]    8
[3656,]    9
[3657,]    8
[3658,]    7
[3659,]    6
[3660,]    5
[3661,]    7
[3662,]    5
[3663,]    4
[3664,]    4
[3665,]    4
[3666,]    4
[3667,]    4
[3668,]    2
[3669,]    3
[3670,]    5
[3671,]    7
[3672,]    7
[3673,]    8
[3674,]    6
[3675,]    6
[3676,]    8
[3677,]    6
[3678,]    7
[3679,]    5
[3680,]    7
[3681,]    7
[3682,]    6
[3683,]    8
[3684,]    8
[3685,]    7
[3686,]    8
[3687,]    6
[3688,]    8
[3689,]    6
[3690,]    8
[3691,]    5
[3692,]    6
[3693,]    5
[3694,]    6
[3695,]    8
[3696,]    8
[3697,]    8
[3698,]    6
[3699,]    8
[3700,]    7
[3701,]    5
[3702,]    5
[3703,]    6
[3704,]    7
[3705,]    3
[3706,]    7
[3707,]    6
[3708,]    2
[3709,]    6
[3710,]    3
[3711,]    7
[3712,]    5
[3713,]    2
[3714,]    7
[3715,]    6
[3716,]    8
[3717,]    7
[3718,]    6
[3719,]    7
[3720,]    6
[3721,]    5
[3722,]    9
[3723,]    6
[3724,]    5
[3725,]    3
[3726,]    6
[3727,]    6
[3728,]    7
[3729,]    4
[3730,]    3
[3731,]    2
[3732,]    5
[3733,]    9
[3734,]    8
[3735,]    6
[3736,]    6
[3737,]    4
[3738,]    3
[3739,]    2
[3740,]    6
[3741,]    8
[3742,]    5
[3743,]    4
[3744,]    5
[3745,]    7
[3746,]    7
[3747,]    5
[3748,]    3
[3749,]    5
[3750,]    7
[3751,]    4
[3752,]    7
[3753,]    5
[3754,]    8
[3755,]    6
[3756,]    7
[3757,]    4
[3758,]    3
[3759,]    4
[3760,]    1
[3761,]    4
[3762,]    6
[3763,]    4
[3764,]    8
[3765,]    9
[3766,]    8
[3767,]    7
[3768,]    7
[3769,]    6
[3770,]    6
[3771,]    5
[3772,]    9
[3773,]    6
[3774,]    5
[3775,]    3
[3776,]    8
[3777,]    4
[3778,]    7
[3779,]    8
[3780,]    7
[3781,]    5
[3782,]    8
[3783,]    4
[3784,]    1
[3785,]    8
[3786,]    6
[3787,]    5
[3788,]    6
[3789,]    6
[3790,]    6
[3791,]    6
[3792,]    7
[3793,]    5
[3794,]    7
[3795,]    6
[3796,]    8
[3797,]    7
[3798,]    1
[3799,]    6
[3800,]    5
[3801,]    6
[3802,]    4
[3803,]    7
[3804,]    4
[3805,]    8
[3806,]    8
[3807,]    7
[3808,]    6
[3809,]    8
[3810,]    3
[3811,]    6
[3812,]    6
[3813,]    7
[3814,]    7
[3815,]    2
[3816,]    7
[3817,]    2
[3818,]    6
[3819,]    3
[3820,]    5
[3821,]    5
[3822,]    7
[3823,]    9
[3824,]    8
[3825,]    5
[3826,]    5
[3827,]    2
[3828,]    5
[3829,]    6
[3830,]    6
[3831,]    6
[3832,]    7
[3833,]    9
[3834,]    5
[3835,]    7
[3836,]    4
[3837,]    5
[3838,]    6
[3839,]    5
[3840,]    6
[3841,]    3
[3842,]    4
[3843,]    4
[3844,]    5
[3845,]    7
[3846,]    8
[3847,]    9
[3848,]    7
[3849,]    7
[3850,]    7
[3851,]    6
[3852,]    6
[3853,]    7
[3854,]    5
[3855,]    6
[3856,]    8
[3857,]    7
[3858,]    9
[3859,]    7
[3860,]    6
[3861,]    3
[3862,]    6
[3863,]    7
[3864,]    7
[3865,]    4
[3866,]    6
[3867,]    7
[3868,]    7
[3869,]    7
[3870,]    7
[3871,]    6
[3872,]    8
[3873,]    6
[3874,]    3
[3875,]    4
[3876,]    6
[3877,]    3
[3878,]    3
[3879,]    6
[3880,]    7
[3881,]    4
[3882,]    6
[3883,]    2
[3884,]    6
[3885,]    1
[3886,]    5
[3887,]    7
[3888,]    6
[3889,]    4
[3890,]    6
[3891,]    4
[3892,]    4
[3893,]    4
[3894,]    9
[3895,]    7
[3896,]    8
[3897,]    4
[3898,]    7
[3899,]    8
[3900,]    7
[3901,]    9
[3902,]    6
[3903,]    4
[3904,]    5
[3905,]    5
[3906,]    4
[3907,]    4
[3908,]    5
[3909,]    7
[3910,]    5
[3911,]    6
[3912,]    6
[3913,]    6
[3914,]    6
[3915,]    3
[3916,]    7
[3917,]    6
[3918,]    6
[3919,]    4
[3920,]    5
[3921,]    5
[3922,]    6
[3923,]    3
[3924,]    4
[3925,]    3
[3926,]    6
[3927,]    8
[3928,]    5
[3929,]    5
[3930,]    5
[3931,]    5
[3932,]    4
[3933,]    8
[3934,]    4
[3935,]    7
[3936,]    8
[3937,]    5
[3938,]    7
[3939,]    6
[3940,]    7
[3941,]    4
[3942,]    5
[3943,]    7
[3944,]    5
[3945,]    6
[3946,]    1
[3947,]    5
[3948,]    6
[3949,]    8
[3950,]    3
[3951,]    6
[3952,]    5
[3953,]    4
[3954,]    3
[3955,]    1
[3956,]    6
[3957,]    5
[3958,]    5
[3959,]    5
[3960,]    2
[3961,]    7
[3962,]    6
[3963,]    5
[3964,]    7
[3965,]    6
[3966,]    7
[3967,]    6
[3968,]    3
[3969,]    6
[3970,]    2
[3971,]    4
[3972,]    8
[3973,]    7
[3974,]    6
[3975,]    6
[3976,]    3
[3977,]    5
[3978,]    7
[3979,]    7
[3980,]    8
[3981,]    7
[3982,]    7
[3983,]    6
[3984,]    5
[3985,]    6
[3986,]    5
[3987,]    4
[3988,]    5
[3989,]    8
[3990,]    4
[3991,]    7
[3992,]    7
[3993,]    8
[3994,]    5
[3995,]    7
[3996,]    8
[3997,]    6
[3998,]    2
[3999,]    7
[4000,]    5
data.frame(obs = ppd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")

Sampling parameters from the prior

Oops, we’ve jumped ahead of ourselves! Best practice is to simulate values from your prior first and check to see if those priors are reasonable.

m1p <-
  brm(data = list(w = 6),                            
      family = binomial(link = "identity"),          
      w | trials(9) ~ 0 + Intercept,                 
      prior(beta(1, 1), class = b, lb = 0, ub = 1),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,          
      sample_prior = "only")
samples_from_prior = as_draws_df(m1p)
samples_from_prior %>% 
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water", title="Prior")

Simulating observations from the prior

We may also want the PRIOR PREDICTIVE DISTRIBUTION which is the expected observiations given our prior.

prior_pd = posterior_predict(m1p)
data.frame(obs = prior_pd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")

Simulating from your priors – prior predictive simulation – is an essential part of modeling. This allows you to see what your choices imply about the data. You’ll be able to diagnose bad choices.

an aside about learning in R

At this point in the course, I’m going to start throwing a lot of code at you. Do I expect you to memorize this code? Of course not.

Do you need to understand every single thing that’s happening in the code? Nope.

But, you’ll learn a lot by taking the time to figure out what’s happening in a code chunk. Class time will frequently include exercises where I ask you to adapt code I’ve shared in the slides to a new dataset or to answer a new problem. When doing so, go back through the old code and figure out what’s going on. Run the code one line at a time. Always observe the output and take some time to look at the object that was created or modified. Here are some functions that will be extremely useful:

str() # what kind of object is this? what is its structure?
dim() # what are the dimensions (rows/columns) of this object
head() # give me the first bit of this object
str(prior_pd)
 int [1:4000, 1] 9 0 0 3 9 8 0 5 4 5 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : NULL
dim(prior_pd)
[1] 4000    1
head(prior_pd)
     [,1]
[1,]    9
[2,]    0
[3,]    0
[4,]    3
[5,]    9
[6,]    8

Continous outcomes

The globe tossing example is cute and easy to work with, but let’s move towards the kinds of variables we more frequently work with. Let’s create a model for some outcome, \(y\) that is continuous.

\[\begin{align*} y_i &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(0, 10) \\ \sigma &\sim \text{Uniform}(0, 5) \end{align*}\]

set.seed(9)
y = rnorm(n = 31, mean = 4, sd = .5)
m2 = brm(
  data = list(y=y),
  family = gaussian,
  y ~ 1,
  prior = c(prior( normal(0,10), class=Intercept),
            prior( uniform(0,5), class=sigma)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.2")
)

An example: weight and height

Using the Howell data (don’t load the rethinking package because it can interfere with brms).

data("Howell1", package = "rethinking")
d <- Howell1
str(d)
'data.frame':   544 obs. of  4 variables:
 $ height: num  152 140 137 157 145 ...
 $ weight: num  47.8 36.5 31.9 53 41.3 ...
 $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
 $ male  : int  1 0 0 1 0 1 0 1 0 1 ...
library(measurements)
d$height <- conv_unit(d$height, from = "cm", to = "feet")
d$weight <- conv_unit(d$weight, from = "kg", to = "lbs")
rethinking::precis(d)
             mean         sd      5.5%    94.5%      histogram
height  4.5362072  0.9055921  2.661042   5.4375      ▁▁▁▁▂▂▇▇▁
weight 78.5079631 32.4502290 20.636856 120.1583 ▁▂▃▂▂▁▁▃▅▇▇▃▂▁
age    29.3443934 20.7468882  1.000000  66.1350      ▇▅▅▃▅▂▂▁▁
male    0.4724265  0.4996986  0.000000   1.0000     ▇▁▁▁▁▁▁▁▁▇
d2 <- d[ d$age >= 18, ]

exercise

Write a mathematical model for the weights in this data set. (Don’t worry about predicting from other variables yet.)

solution

\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

exercise

Simulate from your priors (parameters values and prior predictive values).

solution

Sample from your priors:

m3p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

Sampling parameter estimates for the Intercept.

pairs(m3p)

This is a different (shorter) way to plot your posterior. Good things: it automatically includes the scatterplot so you can see the implications of how these parameters correlate. Bad things: not customizable and not useable when you have a lot of parameters.

Simulate values of weight.

prior_pd = posterior_predict(m3p)
dim(prior_pd)
[1] 4000  352
as.data.frame(prior_pd) %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_histogram() +
  labs(x="Expected observed weights (based on prior)")

Another shorter way:

pp_check(m3p)

Fit the model

m3 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.3"))
posterior_summary(m3)
                Estimate  Est.Error         Q2.5        Q97.5
b_Intercept    99.221909 0.75720356    97.751074   100.737488
sigma          14.291987 0.55198322    13.254363    15.428183
Intercept      99.221909 0.75720356    97.751074   100.737488
lprior         -8.318377 0.05827127    -8.433538    -8.203915
lp__        -1441.294276 1.02695551 -1444.235605 -1440.292326
pairs(m3)
Code
posterior_predict(m3) %>% 
  as.data.frame() %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_density(fill = "grey", color = "white") +
  geom_density( aes(x = weight), data=d2, inherit.aes = F) 
pp_check(m3)

Adding in a linear component

We might assume that height and weight are associated with each other. Indeed, within our sample:

plot(d2$weight ~ d2$height)

exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta h_i \\ \alpha &\sim \text{Normal}(??, ??) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

To update our brms code:

d2$height_c = d2$height - mean(d2$height)
  
m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")
set.seed(9)
samples_from_prior = as_draws_df(m4p)
str(samples_from_prior)
draws_df [4,000 × 9] (S3: draws_df/draws/tbl_df/tbl/data.frame)
 $ b_Intercept: num [1:4000] 167 155.6 149.6 88.7 176.4 ...
 $ b_height_c : num [1:4000] -57.9 -57.2 -18.5 52 -47.3 ...
 $ sigma      : num [1:4000] 14.1 11.8 13.3 2.1 21.1 ...
 $ Intercept  : num [1:4000] 167 155.6 149.6 88.7 176.4 ...
 $ lprior     : num [1:4000] -15.7 -14.7 -12 -15.6 -15.8 ...
 $ lp__       : num [1:4000] -13.8 -12.9 -10.2 -14.9 -14.6 ...
 $ .chain     : int [1:4000] 1 1 1 1 1 1 1 1 1 1 ...
 $ .iteration : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
 $ .draw      : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
Code
d2 %>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_blank() +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_prior[1:50, ],
               alpha=.3) +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))

Describe in words what’s wrong with our priors.

Slope should not be negative. How can we fix this?

Could use a uniform distribution bounded by 0.

m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

exercise

Fit the new weight model to the data.

solution

m4 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4"))
posterior_summary(m4)
               Estimate Est.Error        Q2.5       Q97.5
b_Intercept    99.21508 0.5154139    98.17317   100.20955
b_height_c     24.74661 0.2609169    24.05150    24.99355
sigma          10.36090 0.3898531     9.65234    11.12228
Intercept      99.21508 0.5154139    98.17317   100.20955
lprior        -11.53739 0.0396881   -11.61861   -11.46176
lp__        -1332.15882 1.3936194 -1335.90639 -1330.52002

exercise

Draw lines from the posterior distribution and plot with the data.

solution

Code
set.seed(9)
samples_from_posterior = as_draws_df(m4)
d2%>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_point(size=.5) +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_posterior[1:50, ],
               alpha=.3,
               color="#1c5253") +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))

A side note: a major concern or critique of Bayesian analysis is that the subjectivity of the priors allow for nefarious behavior. “Putting our thumbs on the scale,” so to speak. But priors are quickly overwhelmed by data. Case in point:

m4e = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(-5,5),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4e"))
posterior_summary(m4e)
                Estimate Est.Error         Q2.5       Q97.5
b_Intercept    99.197733 0.5035653    98.214092   100.15932
b_height_c     35.775818 1.8832462    32.177572    39.50070
sigma           9.529429 0.3687178     8.839379    10.27233
Intercept      99.197733 0.5035653    98.214092   100.15932
lprior        -44.172476 3.0738968   -50.456409   -38.49994
lp__        -1334.629214 1.1849503 -1337.645046 -1333.23509

You’ll only really get into trouble with uniform priors that have a boundary, if true population parameter is outside your boundary. A good rule of thumb is to avoid the uniform distribution. We’ll cover other options for priors for \(\sigma\) in future lectures, but as a preview, the exponential distribution works very well for this!